function res = ecogSpiderTSClassify(ecog1,ecog2,c,indexList)
%res = ecogSpiderTSClassify(ecog1,ecog2,c,) Perform a classification on two full ecog datasets
%
% PURPOSE:  Test classifier on dataset. Note that the prediction accuracies 
%           will largely overstimate the accuracies obtainable in a cross  
%           validation!!!! 
%           Only seleceted channels will be included if indexList=[] or omitted.
%           The baseline interval is excluded from the anaylsis if
%           indexList=[] or omitted.
%
% INPUT:
% ecog1:    The first dataset.
% ecog2:    The second dataset. Make sure the same channels are selected as
%           in ecog1.
% c:        OPTIONAL: The penalty paramter. If omitted an estimated is
%           derived from the data.
% indexList: A list of indicex of datapoints to include. This is a list of
%           linear indices as the find command returns it from a matrix of
%           feature weights (e.g. indexList=find(abs(tValues)>criterion);
%           where tValues is a matrix of size channels x timepoints)
%
% OUTPUT:
% res:      A structure holding the results. Fields are
%           loss: The classification error after training 
%           sVs:  Theproportion of support vectors
%           recPos: The recall in the positive class (ecog1)
%           recNeg: The recall in the negative class (ecog2)
%           precPos: The precision in the positive class (ecog1)
%           precNeg: The precision in the negative class (ecog2)
%           nPos:    Then number of examples in the positive class
%           nNeg:    Then number of examples in the negative class
% TODO:     Let user set classification algorithm and parameters 

% 090108 JR wrote it

if exist('svm','file')~=2
    if exist('use_spider','file')~=2
        error('Add the spider toolbox to the matlab path or download it')
    else
        use_spider
    end 
end

if nargin<3
    c=[];
end

if  nargin<4
    indexList=[];
end

% construct input data
% put data in a matrix with dimensions repetitions*features
% on row in features should be ch1(1),...,ch1(end),ch2(1),....ch2(end),...
% MAY REQUIRE TOO MUCH MEMORY. LOOPAROUND TRIALS IF NECESSARY
if ~isempty(indexList)  % a list of inidces to include was passed
    
    data1TrialLength=length(indexList);
    nRepSet1=size(ecog1.data,3);
    nRepSet2=size(ecog2.data,3);
    dataTS=zeros(nRepSet1+nRepSet1,data1TrialLength);
    %the first dataset
    for k=1:nRepSet1
        tmp=ecog1.data(:,:,k);
        dataTS(k,:)=tmp(indexList);
    end
    for k=1:nRepSet2
        tmp=ecog2.data(:,:,k);
        dataTS(nRepSet1+k,:)=tmp(indexList);
    end
    
else    %no indexlist was passed -> exclude baseline, only use selcted channels
    data1TrialLength=prod(size(ecog1.data(ecog1.selectedChannels,ecog1.nBaselineSamp+1:end,1)));
    nRepSet1=size(ecog1.data,3);
    nRepSet2=size(ecog2.data,3);
    dataTS=zeros(nRepSet1+nRepSet1,data1TrialLength);
    %the first dataset
    for k=1:nRepSet1
        tmp=ecog1.data(ecog1.selectedChannels,ecog1.nBaselineSamp+1:end,k)';
        dataTS(k,:)=tmp(:);
    end
    for k=1:nRepSet2
        tmp=ecog2.data(ecog2.selectedChannels,ecog2.nBaselineSamp+1:end,k)';
        dataTS(nRepSet1+k,:)=tmp(:);
    end
end
% Here we have all trials as row vetors in a matrix
objective=ones(size(dataTS,1),1);
objective(nRepSet1+1:end)=-1;   %the labels

% estimate c if necessary
if isempty(c),
    % determine default C according to Joachims
    p=zeros(1,size(dataTS,1));
     for k=1:size(dataTS,1) 
        p(k)=dataTS(k,:)*dataTS(k,:)';
     end
     c=1/mean(p);
end


predictions=zeros(size(objective));
res.losses=zeros(size(objective)); %
sVs=zeros(size(objective));           %collects the number of support vectors
% set classification algorithm
alg = svm;
alg.C = c;
alg.optimizer='libsvm';
d = data(dataTS,objective);
% suppress output
%evalc('[tr resAlg] = train(alg,d)');

[tr resAlg] = train(alg,d);

res.loss = 1-mean(tr.X==tr.Y);
res.sVs=sum(abs(resAlg.alpha)>0)/length(tr.X);
[res.recPos,res.recNeg]=recall(objective,tr.X);
[res.precPos,res.precNeg]=recall(objective,tr.X);
res.nPos = nRepSet1;
res.nNeg = nRepSet2;
res.alg = resAlg;
